function [b , se , p2 , nu] = olsnan_EWC(lhv_full , rhv_full , c , opt)
% This code calculates the standard errors using the equal-weighted cosine
% method. Only deal with number of restrictions = 1 (no F-stat).
% INPUT:
% lhv_full -- Dependent variable, allowed to have nan observations
% rhv_full -- Independent variables, allowed to have nan rows
% c -- whether to include a constant
% opt -- whether use the optimal tuning parameter based on Table II in Lazarus et al. (2018)
%        Default (0): do NOT use estimated autocorrelations and use 0.41
% OUTPUT:
% b -- OLS estimates of coefficients
% se -- standard errors calculated with EWC
% p2 -- two-sided p-values
% nu -- optimal EWC paramter following Lazarus et al. (2018)
%% Initial Clean Up
% Find all non-nan rows
if size(lhv_full , 2) ~= 1
    error('Only univariate LHS allowed')
end

K = size(rhv_full , 2);
a_full = [lhv_full , rhv_full];
a = a_full(sum(~isnan(a_full) , 2) == size(a_full , 2) , :);

lhv = a(: , 1);
rhv = a(: , 2 : end);

T = size(lhv,1);

% If include constant
if c == 1
    rhv_2 = [ones(T , 1) , rhv];
    K = K + 1;
else
    rhv_2 = rhv;
end

% Choice for optimal tuning parameter
switch nargin
   case 3
      opt = 0;
end
%% OLS Estimation
b = olsgmm(lhv , rhv_2 , 0 , -1);
u = lhv - rhv_2 * b;

%% Estimate autocorrelation of regressors and residuals
% We do not do bias adjustment
rho = nan(size(rhv , 2) + 1 , 1);%autocorrelation for each regressor and residual
for j = 1 : size(rhv , 2)
    b_temp = olsgmm(rhv(2 : end , j) , [ones(T - 1 , 1) , rhv(1 : end - 1 , j)] , 0 , -1);
    rho(j) = b_temp(2);
end
b_temp_2 = olsgmm(u(2 : end) , [ones(T - 1 , 1) , u(1 : end - 1)] , 0 , -1);
rho(end , 1) = b_temp_2(2); 

rhos_max = max(rho.^2);%Highest square of autocorrelation

if rhos_max < 0.15
    tune_optimal = 1.62;
elseif rhos_max >= 0.15 && rhos_max < 0.25
    tune_optimal = 1.19;
elseif rhos_max >= 0.25 && rhos_max < 0.35
    tune_optimal = 0.95;
elseif rhos_max >= 0.35 && rhos_max < 0.45
    tune_optimal = 0.78;
elseif rhos_max >= 0.45 && rhos_max < 0.55
    tune_optimal = 0.64;
elseif rhos_max >= 0.55 && rhos_max < 0.65
    tune_optimal = 0.52;    
elseif rhos_max >= 0.65 && rhos_max < 0.75
    tune_optimal = 0.41;    
elseif rhos_max >= 0.75 && rhos_max < 0.85
    tune_optimal = 0.30;
elseif rhos_max >= 0.85
    tune_optimal = 0.18;
end
   
tune_default = 0.41;
%% EWC Parameters
% Optimal choice of free parameter
if opt == 1
    tune = tune_optimal;
elseif opt == 0
    tune = tune_default;
end

nu = floor(tune*T^(2/3));

% Weights 
fvec = pi*(1:nu);
tvec = (0.5:1:T)'/T;
wgt = (sqrt(2)/sqrt(T))*cos(tvec*fvec);  % normalized so that psi'psi = 1

%% Statistics
xu = rhv_2 .* repmat(u , 1 , K);%Time series of X_t*u_t
omega = omega_hat_cos(xu , wgt , nu);
sigma_X_inv = inv((rhv_2'*rhv_2)/T);
varb = 1/T * sigma_X_inv * omega * sigma_X_inv;

se = sign(diag(varb)).*(abs(diag(varb)).^0.5);

t_ratio = b./se;

p2 = nan(K , 1);
for i = 1 : K
    p2(i) = (1 - tcdf(abs(t_ratio(i)) , nu)) * 2;
end
end